Rigorous Multiple-Precision Evaluation of D-Finite Functions in SageMath

We present a new open source implementation in the SageMath computer algebra system of algorithms for the numerical solution of linear ODEs with polynomial coefficients. Our code supports regular singular connection problems and provides rigorous error bounds

NumGfun: a Package for Numerical and Analytic Computation with D-finite Functions

This article describes the implementation in the software package NumGfun of classical algorithms that operate on solutions of linear differential equations or recurrence relations with polynomial coefficients, including what seems to be the first general implementation of the fast high-precision numerical evaluation algorithms of Chudnovsky & Chudnovsky. In some cases, our descriptions contain improvements over existing algorithms. We also provide references to relevant ideas not currently used in NumGfun

Multiple precision evaluation of the Airy Ai function with reduced cancellation

The series expansion at the origin of the Airy function Ai(x) is alternating and hence problematic to evaluate for x > 0 due to cancellation. Based on a method recently proposed by Gawronski, M\"uller, and Reinhard, we exhibit two functions F and G, both with nonnegative Taylor expansions at the origin, such that Ai(x) = G(x)/F(x). The sums are now well-conditioned, but the Taylor coefficients of G turn out to obey an ill-conditioned three-term recurrence. We use the classical Miller algorithm to overcome this issue. We bound all errors and our implementation allows an arbitrary and certified accuracy, that can be used, e.g., for providing correct rounding in arbitrary precision

Rounding Error Analysis of Linear Recurrences Using Generating Series

We develop a toolbox for the error analysis of linear recurrences with constant or polynomial coefficients, based on generating series, Cauchy's method of majorants, and simple results from analytic combinatorics. We illustrate the power of the approach by several nontrivial application examples. Among these examples are a new worst-case analysis of an algorithm for computing Bernoulli numbers, and a new algorithm for evaluating differentially finite functions in interval arithmetic while avoiding interval blow-up

The expected signature of Brownian motion stopped on the boundary of a circle has finite radius of convergence

The expected signature is an analogue of the Laplace transform for rough paths. Chevyrev and Lyons showed that, under certain moment conditions, the expected signature determines the laws of signatures. Lyons and Ni posed the question of whether the expected signature of Brownian motion up to the exit time of a domain satisfies Chevyrev and Lyons' moment condition. We provide the first example where the answer is negative

Computing roadmaps in smooth real algebraic sets

International audienceLet (f1, . . . , fs) be polynomials in Q[X 1 , . . . , Xn ] of degree bounded by D that generate a radical equidimensional ideal of dimension d and let V ⊂ C^n be the locus of their complex zero set which is supposed to be smooth. A roadmap in V ∩ R^n is a real algebraic curve contained in V ∩ Rn which has a non-empty and connected intersection with each connected component of V ∩ R^n . The classical strategy to compute roadmaps is due to J. Canny and leads to algorithms having a complexity within D^O(n^2) arithmetic operations in Q. This strategy is based on computing a polar variety of dimension 1 and a recursion on the studied variety intersected with fibers taken above a critical value of a projection. Thus, it requires computations with real algebraic numbers and introduces singularities at each recursive call. Thus, no efficient implementation of roadmap algorithms have been obtained until now. Our aim is to provide an efficient implementation of the roadmap algorithm. We show how to slightly modify this strategy in order to avoid the use of real algebraic numbers and to deal with smooth algebraic sets at each recursive call in the case where the input variety is smooth. Our complexity is h^d D^O(n) operations in Q where h bounds the number of recursive call in our algorithm. This quantity is related to the geometry of V ∩ R^n and is bounded by D^O(n), thus in worst cases our algorithm has a complexity within D^O(n^2) arithmetic operations. We report on some experiments done with a preliminary implementation of our algorithm

Comparison between binary and decimal floating-point numbers

International audienceWe introduce an algorithm to compare a binary floating-point (FP) number and a decimal FP number, assuming the "binary encoding" of the decimal formats is used, and with a special emphasis on the basic interchange formats specified by the IEEE 754-2008 standard for FP arithmetic. It is a two-step algorithm: a first pass, based on the exponents only, quickly eliminates most cases, then, when the first pass does not suffice, a more accurate second pass is performed. We provide an implementation of several variants of our algorithm, and compare them